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A general thermal ablation model for silicates is proposed. The model includes the mass losses through 
the balance between evaporation and condensation, and through the moving molten layer driven by surface 
shear force and pressure gradient. This model can be applied in the ablation simulation of the meteoroid 
and the glassy ablator for spacecraft Thermal Protection Systems. Time-dependent axi-symmetric 
computations are performed by coupling the fluid dynamics code, Data-Parallel Line Relaxation program, 
with the material response code, Two-dimensional Implicit Thermal Ablation simulation program, to 
predict the mass lost rates and shape change. The predicted mass loss rates will be compared with available 
data for model validation, and parametric studies will also be performed for meteoroid earth entry 


conditions. 


Nomenclature 


= normalized ablation rate 

= specific heat, J/kg-K 

= heat transfer coefficient, J/m?-s-K 
= mass transfer coefficient, kg/m?-s 
= diffusion coefficient, m7/s 

= total internal energy, J/m* 

= enthalpy, J/kg 


_ Pvby—Pchc 


, J/k 
Pv-Pc 





“ Aerospace Engineer, Thermal Protection Materials Branch, MS 234-1. 


American Institute of Aeronautics and Astronautics 


H, = recovery enthalpy, J/kg 


k = thermal conductivity of solid, W/m-K 

M = molecular weight, kg/kmol 

m = mass flux, kg/m?-s 

Dp = pressure, Pa 

p = saturated vapor pressure, Pa 

q = heat flux, W/m? 

Yeond = in-depth conductive heat flux at surface, W/m? 


Uconv = convective heat flux at surface, W/m? 

Irw = radiative heat flux at surface, W/m? 

Ory = translation-vibration energy exchange rate, W/m? 
r = radius, cm 

fe = corner radius of the model, cm 

Tn = nose radius of the model, cm 

R = universal gas constant, 8.3143 J/ mol-K 

T = temperature, K 

AT rex = maximum allowed change on surface temperature, K 
T,, = environment temperature, K 

u = flow velocity, m/s 

v = speed of moving grid, m/s 

Vv = thermal velocity, m/s 

x,y = Cartesian coordinate system, m 

w = source term of gas-surface interactions, kmol/m?-s 
a = absorptance 

Y = evaporation/condensation sticking coefficient 

0 = effective thickness of moving liquid layer, cm 
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é = emissivity 


p = total density, kg/m* 

X = blowing reduction coefficient 
uu = viscosity, Pa-s 

o = Stefan-Boltzmann constant, W/(m?-K*) 
T = shear force, Pa 

subscripts 

1 = unblown 

c = char 

g = pyrolysis gas 

ij = direction components 

l = molten layer 

n = surface species 

S = gas species 

Vv = vapor or vibrational energy 
w = wall 


I. Introduction 


Chondrites are the most common type of meteoroid that falls to Earth. Most chondrites are rich in silicate materials. 
To determine the physical and chemical nature of the meteoroid, one must understand the interaction between the 
meteoroid and the atmosphere. In the classic meteor ablation model [1], ablation begins when the surface of the 
meteoroid reaches the boiling temperature, and then the temperature is assumed to remain constant. Classic meteor 
ablation theory assumes that the mass loss is proportional to the kinetic energy imparted to the meteoroid. Meteoroids 
may lose mass not only though evaporation, but also through spraying of the molten layer on the surface. To simulate 
this effect, the classic ablation model includes a mass loss term which is an inverse tangent function of the difference 
between surface temperature and melting point. The other commonly used mass loss model for meteoroid is the mass 


loss rate computed using Langmuir evaporation assuming the condensation rate is negligibly small. [2] 
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In this work, we propose a general thermal ablation model for silicate materials. The proposed ablation model 
includes the mass losses through the balance between evaporation and condensation, and through the moving of molten 
layer driven by surface shear force and pressure gradient. This proposed ablation model is implemented in the Two- 
dimensional Implicit Thermal AblatioN simulation program, TITAN [3], which was developed for the computation 
of Thermal Protection System materials under hypersonic entry conditions. The ablation characteristics of silicates 
are the same as those of silica and metal, which are the commonly used materials for the Thermal Protection Systems 
for NASA’s space missions. Hence, this ablation model can be applied for both the meteoroid and the glassy or 
metallic components of spacecraft Thermal Protection Systems. The time dependent axi-symmetric simulations for 
mass loss and shape change predictions are performed by coupling the fluid dynamics code, Data-Parallel Line 
Relaxation (DPLR) [4], and the material response code, TITAN. The predictions will be compared with available data 
for model validation. [5] Parametric studies under meteoroid entry conditions will be performed to gain further 


understanding on the ablation response of silicate materials to guide the direction of future work. 


II. Thermal Ablation Model 


The mechanism of thermal ablation for silicate materials includes the motion of molten layer driven by surface 
shear force and surface pressure gradient, m , and the balance between evaporation and condensation at the surface, 
My. 

When subjected to aerothermodynamics heating and shear forces, the molten material will form a flowing liquid 
layer over the solid surface. The equations governing the liquid layer are nearly the same as for the case of an 
incompressible gas boundary layer, the difference being a variation of viscosity with temperature. For an axi- 
symmetric geometry, if the inertia force in momentum equations can be neglected, the momentum and continuity 


equations of molten layer can be integrated to obtain the mass loss rate due to the motion of liquid layer [6]: 


rin) =" (1,0 — 2x 53)) (1) 


r ax “Uw dx 


The viscosity of molten layer usually increases exponentially with a decrease in temperature i.e. 
a 
wser™ (2) 
Here a and b are constants. The “effective thickness” of a moving liquid layer, 0, is defined as the distance at which 


the viscosity increases by one exponential factor and it’s assumed that the thickness of the moving liquid layer is thin 


in comparison with the thickness of thermal layer. 
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The evaporation rate of a surface species can be computed from its saturated vapor pressure, P’, and its 
condensation rate from the gas phase can be computed if the species partial surface pressure, P, , is known [7] : 
tity = Dn Mn (PE — Pa) (3) 
Here v,, is the thermal speed of surface species n, and y is the evaporation/condensation (or sticking) coefficient. 
If the ambient gas pressure is very low, the condensation rate is thus negligible. Then the computation of m, can 
be simplified to 
Ty = Dn LYM, Pe (4) 
The saturated vapor pressure can be computed using a chemical equilibrium code, such as MAGMA [8]. The 
MAGMA code was developed specifically for the composition of meteoroids or earth crust. Forsterite, Mg2SiO«(s), 
is one of the silicates which we are interested in this study. n MAGMA, Mg>SiO,(1) is assumed to undergo pseudo- 
reaction to form pseudo-species, MgO(1) and SiO2(1) which then becomes vapor, MgO(g) and SiO2(g), as the 
temperature increases. The computed saturated vapor pressures as function of surface temperature for MgO(1) and 
Si02(1) are presented in Fig. La. 
If a continuous flow is formed over the surface, the species partial pressure at the surface has to be calculated using 
a finite-rate reacting flow solver. The gas phase chemical species and their chemical reactions for ablation products 
have to be considered in the reacting flow simulation [9]. A finite-rate reacting flow simulation with all the chemical 
species of ablations products is complicated and extremely CPU-time consuming. However, if the chemical 
equilibrium condition can be established at the surface, the net mass loss rate due to evaporation/condensation process 
can be pre-calculated, i.e. the B’ tables, using a chemical equilibrium code, such as MAT [10]. The computed B¢ of 
Mg2SiO, as a function of surface temperature and pressure are presented in Fig. 1b. The mass loss due to 
evaporation/condensation for a chemical equilibrium surface is written as 
My = CyBe (5) 
The total mass loss rate due to thermal ablation, ™,, is thus expressed as 


Me = My +m, (6) 
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Figure la: Saturated vapor pressures for SiO2 and MgO. 
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Figure 1b: The B’ tables at various surface pressures for Forsterite. 
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III. Computational Methods 


In this section we describe the relevant equations of the TITAN and DPLR codes and the coupling procedure 

employed to obtain the solutions presented in this work. 
A. Material Ablation and Thermal Response Code 

The governing equations in TITAN, which include energy conservation, a three-component decomposition model, 
and the surface energy balance, are solved with a moving grid system to simulate the shape change due to surface 
recession [3]. The equations are converted from a Cartesian coordinate system (x,y) into a general body-fitted 
coordinate system (é, 7) and then discretized using a finite-volume method. A moving grid technique was 
implemented to solve simulations with surface recession. 

The internal energy balance is a transient thermal conduction equation with additional pyrolysis terms: 

Ply 5 = V- (KV) — (hg —h)V- ting — tity Why + peyv VT (7) 

The individual terms in Equation (1) may be interpreted as follows: the rate of storage of sensible energy, the net rate 
of thermal conductive heat flux, the pyrolysis energy-consumption rate, the net rate of energy convected by pyrolysis 
gas, and the convection rate of sensible energy due to coordinate system movement. There is no pyrolysis gas flow 
inside the silicate materials, m, = 0, and thus the second and the third terms in Equation (7) are equal to zero. 

Conditions at the ablating surface are determined by the aerothermal environment and by chemical interactions 
between the boundary layer gas, the pyrolysis gas, the ablation products, and the chemical constituents of the surface 
material. TITAN employs a convective transfer-coefficient form of the surface energy balance: 

Cu (Hy — Ay) + tithe + tighg — (me + Titg hy + Mw 4rw — FEw(Tw — TS) — Geona = 0 (8) 

Here p.ueCy is the convective heat transfer coefficient, H, is the recovery enthalpy, and all other quantities are defined 
at the ablating surface. The first term in Equation (8) is the convective heat flux, the second through fourth terms 
represent the chemical energy released (or absorbed) by ablation, the fifth and sixth terms are radiation absorption and 
emission, respectively, and the final term is the rate of heat conduction into the TPS. This equation has been simplified 
from the general form by assuming equal diffusion coefficients within the boundary layer and equal Stanton numbers 
for heat and mass transfer. These are standard assumptions for most entry environments and TPS materials. If these 
assumptions are not applicable, a more general form of the surface energy balance may be used. Flow-field radiation 


is reflected or absorbed at the surface, but not transmitted. 
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A blowing correction accounts for the reduction in heat transfer coefficient due to the injection of gases from 
pyrolysis and surface ablation into the boundary layer. The blowing correction equation used by TITAN is: 
(mc+Mg) 


Cc In(1+24B’ 
-H a where B' = 
CH1 2AB Cu 


(9) 


Here / is the blowing reduction parameter, and Cy/Cy1 is the ratio of the blown (ablating) to the unblown (non- 
ablating) heat transfer coefficients. For laminar flow, / is 0.5 or higher, depending on the geometry and the ratio of 
molecular weights of the injected and boundary-layer-edge gas. For transitional or turbulent flow, smaller values of / 
are used. B’ is the non-dimensional mass blowing rate at the surface. Unless noted otherwise, a blowing reduction 
parameter of 0.5 is used for the calculations presented herein. By assuming that the momentum boundary layer 
thickness is the same as the thermal boundary layer thickness, the surface shear force for an ablating surface is 


approximated by the following equation: 


Sr (10) 


Tw1 = C1 


B. Fluid Dynamics Code 

The Navier-Stokes solver, DPLR, is used to estimate the hypersonic aero-heating distribution over a blunt body 
[4]. The governing equations may be characterized as representing a flow-field in thermal and chemical non- 
equilibrium. The DPLR code solves the time-dependent conservation equations of mass, momentum, and energy 


within the flow-field. The species mass conservation equation is given by: 


dps . Oa a 
FE Bag (Pty) = ~ ayy (Peas) + Ms ay 


The total momentum conservation is written as: 


OT; 


a a 
5 (Pui) + ax; (pujuj) = — Gay (12) 
The vibrational and total energy equations are written as: 
dE, , a _ a 
“att ax (Foy) = — 5 (aj) + Orv (13) 
dE, a a a a 
ra ax, (( + p)u) ae (4; + Wj) - ax (Mitis) = Ls=1 53, Dsjhs (14) 


In DPLR, the Gauss-Seidel line relaxation method is modified to enable fast convergence for viscous flow with 


high parallel efficiency on massively parallel computers. 
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C. Fluid-Solid Coupling 

Coupling between the material response and flow codes is required for many multidimensional ablation 
simulations, because the magnitude and distribution of the surface heat flux, pressure, and shear force are very 
sensitive to shape change. The integration between TITAN and DPLR is based on a loosely coupled approach 
following the flowchart presented in Fig. 2 [11]. In this approach, all CFD solutions are computed for an unblown 
surface but use the ablated shape. Chemical equilibrium surface thermochemical interactions and blowing effects are 
incorporated in the material-response code by use of ablation tables, the surface energy balance with heat transfer 
coefficient, and the blowing reduction parameter, as described by Equations (8) and (9). 

The initial flow field and its associated surface heat flux, pressure, and shear are first computed by DPLR. For 
each surface point, heat transfer coefficient, Cy;, shear force, T,,;, and pressure, p, are calculated. Then shear force, 
pressure, and heat transfer coefficient are passed as boundary conditions to TITAN. A time-accurate ablation and 
thermal response computation is performed by TITAN. When the maximum surface temperature or recession meets a 
pre-specified limit, TITAN stops its computation and outputs the location of the ablated surface. A new CED grid is 
then generated based on this ablated surface and estimated shock location, and a new steady-state flow solution is then 
calculated by DPLR. The pressure, shear force, and heat transfer coefficients are calculated for this new solution. The 
updated values for pressure, shear force, and heat transfer coefficients are input to TITAN for another run of time- 
accurate ablation and in-depth thermal response. This procedure is repeated until TITAN reaches the specified final 
time. In each CFD run, the outer boundary of grid is aligned with the shock using an internal subroutine in DPLR, and 
the cell Reynolds number equal to 1 is enforced. The TITAN-DPLR coupling procedure is written in a UNIX/Linux 


script file. The test cases presented in this paper are computed on a Linux cluster system. 


American Institute of Aeronautics and Astronautics 





TITAN 


(time marching from 
t* to t*+ At,) 









AT,, @t*+Aty < AT ma 







No DPLR 


(steady-state) 












Update BC’s at] | = 
(t+ At 


Figure 2. DPLR/TITAN shape change coupling methodology. 


IV. Results 


The computations presented in this section include two parts. In the first part, ablation model validation will be 
performed. The ablation characteristics of silica (SiOz) are the same as those of silicate materials, and the thermal 
properties of quartz are well studied. Thus, quartz is one of the best candidate materials for our ablation model 
validation. The ablation data for quartz published in the work of Adams et al. [5] will be used to evaluate the accuracy 
of our ablation models. The preliminary results performed by coupling DPLR and TITAN for a 10-deg half angle 
Opaque quartz sphere-cone are presented in Figures 3. The history of surface contours is shown in Figure 3a, and the 
temperature contours at 25 seconds are presented in Figure 3b. Comparisons of stagnation point recession vs time 
between computation and data for the opaque and clear quartz models will be presented in the final paper. In the 
second part, the ablation simulations and parametric studies for Forsterite (Mg2SiO«) and typical chondritic material 
will be discussed. The flow environments used in the second part of computation are typical for the meteor earth entry. 


The maximum surface heat flux is in the order of MJ/cm?, and the maximum surface pressure is around hundreds bars. 
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Figure 3a: Shape change as function of time 
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Figure 3b: Temperature contours at time equal to 25 seconds 
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